Defect Structures in the Growth Kinetics of the Swift-Hohenberg Model 
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The growth of striped order resulting from a quench of the two-dimensional Swift-Hohenberg 
model is studied in the regime of a small control parameter and quenches to zero temperature. We 
introduce an algorithm for finding and identifying the disordering defects (dislocations, disclinations 
and grain boundaries) at a given time. We can track their trajectories separately. We find that the 
coarsening of the defects and lowering of the effective free energy in the system are governed by a 
growth law L(t) ~ t x with an exponent x near 1/3. We obtain scaling for the correlations of the 
nematic order parameter with the same growth law. The scaling for the order parameter structure 
factor is governed, as found by others, by a growth law with an exponent smaller than x and near 
to 1/4. By comparing two systems with different sizes, we clarify the finite size effect. We find that 
the system has a very low density of disclinations compared to that for dislocations and fraction of 
points in grain boundaries. We also measure the speed distributions of the defects at different times 
and find that they all have power-law tails and the average speed decreases as a power law. 

PACS numbers: 05.70.Ln, 64.60. Cn, 64.75.+g, 98.80.Cq 
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I. INTRODUCTION 

What are the defects which control the long-time or- 
dering of systems growing a striped pattern? This ques- 
tion arises in a variety of physical contexts [1]. Here we 
are motivated by the recent experiments [2,3] investigat- 
ing the ordering of a two dimensional diblock copolymer 
system. The system studied offers a physical realization 
of the ordering in an isotropic two-dimensional smectic 
material. In these experiments they found that the late- 
time ordering satisfies scaling with a growth law L « t x 
with x = 1/4 and the final stages of ordering are governed 
by the annihilation of sets of disclination quadrapoles. In 
this paper we address the question: Is the ordering in this 
physical system described by the Swift-Hohenberg (SH) 
model [4] , the simplest model one can construct to govern 
the ordering in stripe forming systems? 

We investigate the growth kinetics of the Swift- 
Hohenberg model for a small control parameter (e = 0.1) 
in two dimensions and quenches to zero temperature. It 
is this regime which appears most likely to correspond to 
the experimental situation. In large e regime the system 
evolves to a glassy state. We focus primarily on the de- 
fect structures generated in the ordering of the system. 
In the most naive picture of this ordering process one can 
think in terms of an initial local layering, as in a smec- 
tic, in some direction. This ordering can be disrupted by 
point defects: dislocations and disclinations. This sug- 
gests a coarsening picture with annihilating point defects 
similar to the case of the XY model [5] and a growth law 
with exponent x = 1/2. This simple picture is not seen 
in simulations. We find, in agreement with the numerical 
results of Hou et al. [6] and Boyer and Vinals [7], that 
the defect structures for the SH model arc dominated by 
grain boundaries which persist for long times. Unlike the 
case of an XY model, the ordering is not dominated by 



annihilation of isolated point defects. These are observed 
but are not the dominant structures. 

We find numerically, at late times after finite size ef- 
fects enter, that the system becomes anisotropic, and the 
grain boundaries shrink. In this case one sees a cross-over 
to an effective growth exponent x = 1/2. 

We give below a detailed numerical study of the statis- 
tical properties of the defects disrupting striped pattern 
formation in the SH model. In order to carefully discuss 
the defects we need a reliable filter for finding them. We 
present an algorithm which effectively locates defects and 
grain boundaries for any control parameter e. We can 
distinguish between grain boundaries and other defects, 
and track their trajectories separately. We compare this 
method to the other approaches used in earlier work in 
appendix A. 

There are a number of ways of characterizing the de- 
gree of ordering in these systems: (i) Counting the num- 
ber and size of defects and their evolution with time, 
(ii) Monitoring the lowering of the average effective driv- 
ing free energy as a function of time, (iii) Evaluation of 
the nematic order parameter correlation function and its 
associated scaling behavior, (iv) Evaluation of the or- 
der parameter structure factor and its associated scaling 
behavior. We find that (i), (ii), and (iii) can all be char- 
acterized by a single growth law with the exponent near 
1/3, while the order parameter scaling, as found by oth- 
ers, is characterized by a growth law with the exponent 
near 1/5. 



II. SWIFT-HOHENBERG MODEL 



The Swift-Hohenberg model for a scalar order param- 
eter, ip, is specified by the equation of motion 
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where e is a positive control parameter, q n is the magni- 
tude of an ordering wavenumber and £ is the Gaussian 
noise satisfying (C(r, t)((r' , t')) = 2T5{v - r')S(t - t 1 ), 
where the noise strength T is proportional to the final 
temperature governing the system after a quench. We 
will focus here on quenches to zero temperature where 
we can set T and the noise ( to zero. We are interested 
in the growth kinetics problem where we prepare this sys- 
tem initially in a completely disordered state. We then 
allow the system to evolve forward in time to form a 
striped pattern. For example one could choose 



(V(x,t )V(y,io)> = *§5(x-y) 



(2) 



where is a constant. However the precise form of the 
initial conditions is not important [5]. 

This model can be formulated as a Langevin equation 
driven by an effective Hamiltonian: 



n E = j '^{-|v 2 + ^ 4 + ^[(% 2 + v 2 )v] : 



If we introduce 



E(t) = (H E )t 



(3) 



(4) 



where the average is over an ensemble of initial condi- 
tions, then E(t) is lowered as the system orders in a 
striped pattern with wavenumber q^. 

Eventually the system approaches an ordered state de- 
scribed approximately by the single-mode approximation 
[8] where, assuming layering along the z-direction, 



■00 = Aq cos q z . 



(5) 



If we put this ansatz into Eq.(3), assume that the system 
is an integral number of wavelengths in the z-direction, 
and minimize with respect to the amplitude Aq, we ob- 
tain the results, 
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where S is the area of the system. Pomeau and Man- 
neville [8] have shown that this is a very good approxi- 
mation for the "ground" state even for moderately large 
values of e. In the growth kinetics context the approach 
to equilibrium is monitored by 



AE(t)=E(t)-E eq <xL E \t) , 



(9) 



and 



^ 2 (t)^(^l)-(^) t ^L-\t) , 



(10) 



where oc L^(t) [9]. 

Another measure of the ordering in the system is given 
by considering the director field 



n(x) = 



W(x) 



|VV(x)| ' 

and the associated nematic order parameter 

1 



(11) 



Qa/3 — Qo 



n a n/3 - - 5 af 3 



(12) 



In two dimensions, however, all of the information in this 
order parameter is contained in the quantity cos 29 where 
h = (cos 9, sin 0). It is easy to show, for example, that 

C nn (x,y,t) =2<TrQ(x,t)Q(y,t)) t 

= <cos[( V (x,t)- V (y,t)]) t . (13) 

where 



If we define 



then 



(p(x, t) = 20(x, t) . 



B x = n 2 x - hi 
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C nn (x,y,t) = (B(x,t)-B(y,t)) t 



(14) 



(15) 



(16) 



(17) 



The nematic order parameter correlation function, 
C„„, was shown by Christensen and Bray [10] to obey 
scaling in the conventional form 



C nn (r,t)=F(r/L n (t)) , 



(18) 



where r = x y. Elder, Vinals and Grant [11] showed 
that the scaling of the order parameter structure factor 

S(k,t) = (\Mt)\ 2 ) = Ls^F^k - q Q )L s (t)) , (19) 

differs from that observed in ordering system without 
stripes: S(k,t) = L 2 (t)F 2 (kL(t)). 
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III. REVIEW OF PREVIOUS WORK 

The early work on this problem focused on establishing 
the final equilibrium state reached after a quench. This is 
a two dimensional system and by forming stripes one has 
a broken continuous symmetry. The behavior of the sys- 
tem at non-zero temperatures, as for the two dimensional 
X-Y model, requires, as pointed out by Toner and Nelson 
[12] a treatment of both long wavelength fluctuations in 
the layers and free defects. Above a Kosterlitz-Thouless 
type transition one has an isotropic phase while below 
this transition one has a phase with persistent orienta- 
tional order. 

In an early paper, Elder, Vihals, and Grant [11] car- 
ried out a numerical analysis leading to the scaling solu- 
tion given by Eq.(19). Working with fixed e = 0.25 they 
looked at the system's ordering as a function of noise 
strength T. They found a qualitative difference between 
low noise and high noise. For the large noise case they 
found a rapid (exponential) relaxation to the asymptotic 
stationary state and a power-law approach for the lower 
noise case. Their results are in agreement with the pic- 
ture due to Toner and Nelson that one has a transition 
to an isotropic state for large enough noise. There is no 
real ordering in the isotropic state and this is why there 
is exponential decay to the equilibrium state. In the or- 
dered state one has scaling and a power-law growth law 
which, for small noise, they found to have an exponent 
x s = 1/4. They found a smaller exponent x s = 1/5 at 
low temperatures, but they had less statistics and there 
appeared to be " difficulty removing defects" . They ar- 
gued for a late time cross over to the expected x = 1/2 
but they did not see this. 

Cross and Meiron [13] also studied the SH model nu- 
merically in the absence of noise. They found a x s = 1/4 
for e = 0.25. The dynamics appear to freeze for higher e. 
They looked at the defect structure but in a qualitative 
way noting the existence of domain walls rather than a 
set of isolated point defects. The theoretical discussion 
in their paper is based on the phase-field approximation 



dj)_ 
~dt 



(D||Vjj + D±Vi 



(20) 



which from the most naive point of view suggests a 
growth law with exponent x = 1/2. They discuss some 
selection mechanisms which could lead D\\ and to 
adjust themselves to zero and reduce x to 1/3 or 1/4. 
They concluded that they did "not have a good theoret- 
ical understanding of these results" and suggested that 
the defects in the problems should be treated explicitly. 

Hou, Sasa, and Goldenfcld (HSG) [6] confirmed previ- 
ous numerical results which showed for e = 0.25, x s = 1/5 
with zero noise and x s — 1/4 with nonzero noise as ob- 
tained from the structure factor scaling. They went fur- 
ther and used a simple method to identify domain walls 



and measure their lengths (more about this below) . They 
measured excess energy, AE(t), and the domain wall 
length and found that they show the same scaling ex- 
ponents 1/4 at zero noise and 0.3 at non-zero noise. The 
energy does go to the lowest order in e value of — e 2 /6 
in the noiseless limit. They find "defects are indeed the 
driving force behind the coarsening process due to its 
dominant contribution to the excess energy." They sug- 
gest that the phase field approach gives the wrong ex- 
ponent because it does not include the effects of defects. 
For larger e (=0.75) they found much slower logarithmic 
growth. The system seems to become glassy. 

Christensen and Bray [10] also carried out numerical 
work on the SH model for e = 0.25 and found x s = 1/5 
for zero noise and x s — 1/4 for nonzero noise. From 
scaling of the director correlation function they find ex- 
ponents are 0.25 and 0.30 for zero and nonzero noise. 
They suggest that there is a cross over to x = 1/2 at 
very long times. The theory they developed does not 
include defects. 

Boyer and Vihals [7] point out "Near the bifurcation 
threshold, the evolution of disordered configurations is 
dominated by grain boundary motion through a back- 
ground of largely immobile curved stripes" . They find 
for small e an exponent x = 1/3 which they interpret as 
arising from a law of grain boundary motion [14] . Else- 
where [15] they also point out for larger values of e the 
dynamics cross over to a frozen state with quenches to 
zero temperature. This glassy behavior is associated with 
grain boundary pinning. 



IV. NUMERICAL RESULTS FOR SH MODEL 

We present here our numerical results for the SH equa- 
tion. We follow the numerical prescriptions of Bray and 
Christensen [10]. We use the finite difference scheme on 
two dimensional lattice of sizes 256 x 256 and 512 x 512 
with periodic boundary conditions. We set e = 0.1, 
Ar = tt/4 and At = 0.03. We replace d t ip(r,t) by 
« +1 " Wj) 1^ and V 2 VM) by 
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(Ax) 2 



NN 



NNN 



(21) 



where NN and NNN mean the nearest neighbors and 
next-nearest neighbors respectively. By choosing the 
proper scale of time, space and the order parameter, we 
can set qo = 1. The systems have eight grid points per 
wavelength. We used uniformly distributed random ini- 
tial conditions. 

For the smaller 256 x 256 systems we were able to fol- 
low the ordering process to very late stages. Some of the 
independent trials proceed to a final state where we have 
a set of well aligned layers. 
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In Fig. 1 we plot AE(t) and Aip 2 (t) for an ensem- 
ble of runs on a 256 x 256 lattice. We note that there 
are two regimes where Le defined by Eq.(9) is described 
by different exponents. For t s < t < t c (£ s w 300 and 
£ c w 9000) we find xe ~ 0.3, while for £ > t c we find 
xe ~ 0.5. The cross over at £ s=s f c appears to be due to 
the finite size effects, as we discuss below. For £ > £ c the 
system is effectively anisotropic and we find an effective 
exponent xe near to 1/2. 
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FIG. 1. AE(t) and Aip 2 (t) for a 256 X 256 system. Straight 
lines are used to fit different parts. Averaged over 40 trials. 



In Fig. 2 we plot AE(t) and Ai/> 2 (£) for a 512 x 512 
system. In this case we see that £ c has been extended to 
much larger values and we have not been able to follow 
the ordering process to completion. Our fits to A_E(£) 
and Aip 2 it) in the regime £ s < £ < £ c again gives, to 
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FIG. 2. AE(t) and Ai/> 2 (t) for a 512 X 512 system. The data 
for t < t s are not shown. The straight lines are used to guide eyes. 
Averaged over 57 different trials. 

To probe directly the stripes' increasingly orientational 
order, we measure the nematic order parameter correla- 
tion function C nn (r, t) in the 512 x 512 system. The re- 
sults, averaged over 57 runs, are shown in Fig. 3. We ob- 
tain scaling with a correlation length obeying the growth 
law L n oc £ - 36 . We can estimate the time £ c when the 
cross-over begins in this larger system as follows. The 
system becomes anisotropic and one expects cross-over 
when the correlation length L n grows to be some sub- 
stantial fraction of a lateral dimension of the system. In 
terms of ratios we can write 



L»(£ c (512)) _ 512 
L„(£ c (256)) ~ 256 



£ c (512) 
£ c (256) 



1/3 



(22) 



higher accuracy, xe — x^ = 1/3. 



In the 256 x 256 system £ c (256) ~ 9000, so we obtain 
£ c (512) ~ 60000. Notice that in Fig. 2 the effective ex- 
ponent xe begins to increase at the time 50000 ~ 70000, 
which is consistent with our estimate. 
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FIG. 3. Time evolution of the correlation function C nn (r,t) in 
512x512 SH system illustrated with times 6x 10 3 , 1.2x 10 4 , 1.8xl0 4 
increasing from left to right. We extract the time evolution of 
the correlation length L(t) by monitoring the r a (t) for which 
C nn (r a (t)) = a, where we choose a = {0.3,0.4,0.5,0.6,0.7}. The 
scaling exponent x n is extracted from the log-log plot insert of r a (t) 
v.s. t by fitting it with a straight line. Averaged over 57 trials. 
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FIG. 4. The structure factor S(k,t) = (|i/>(k, t)| 2 ) in the 
512 X 512 system. The log-log plot of S(qo,t) v.s. t can be fit 
to t x with x = 0.24. The scaling collapse of the structure factor 
was obtained with x = 0.24 as the scaling exponent. Averaged over 
57 independent trials. 

In Fig. 4 we plot the structure factor S(k, t) and show 
that scaling holds in the form given by Eq.(19) with a 
growth law characterized by an exponent x s = 0.24 as 
shown in the insert. Our results here agree with those 
found previously that the exponent governing the growth 
law for the structure factor is significantly smaller than 
that governing the nematic order parameter. 



V. DEFECT STRUCTURES AND DYNAMICS 

In Fig. 5 we show a typical configuration for the Swift- 
Hohenberg model for a quench to zero temperature after 
a time 12000 for a 512 x 512 system. Notice the rather 
complicated structure which includes dislocations, discli- 
nations and grain (domain) boundaries. Our main focus 
in this paper is to study the statistics of these defects. In 
appendix A we discuss an algorithm for picking out the 
defects and tracking their motions. 




FIG. 5. A typical configuration for the SH model for a quench to 
zero temperature after t = 12000 in a 512 X 512 system. The black 
points correspond to V( x ) > 0, and the white points to V( x ) < 0. 

If we look at Fig. 5 we see that it shows a compli- 
cated situation with a variety of different defect struc- 
tures which one can identify by eye at the length scale 
of several layer spacings. At a more fundamental level 
we need a way of identifying which points in space, at 
the level of each site on the numerical grid, are part of 
a defect. At the shortest length scale in the problem the 
order parameter is Q a p defined by Eq.(12). For this two- 
dimensional system this can be replaced by the vector 
order parameter B defined by Eqs.(15) and (16). The 
assumption is that all of the defects in the system can be 
built up from the ±^ disclinations in the director field n 
which translate into vortices with charge ±1 for the field 
B. We identify these defects by looking for the cores of 
the vortices. We can find the cores of the defects by look- 
ing from those sites where B is changing rapidly. We can 
define 



(23) 



and identify defect points as those sites where A is larger 
than some value. Notice that A can also be written in 
the form 
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,4 = 4^(V Q n (3 ) 2 = (V Q ^) 2 . (24) 

a, /3 

The precise numerical determination of A is discussed in 
appendix A. Notice that A is proportional to the gradient 
energy for an isotropic nematic. 

In analyzing their experimental data Harrison, et al. 
[2,3] found a set of fundamental disclinations and from 
these built up dislocations as bound disclinations with 
opposite charge. They used this procedure to identify a 
large dislocation density. Most of the fundamental discli- 
nations went into forming these dislocations since in the 
end the ratio of dislocations to the remaining disclina- 
tions was about ten to one. In our case the situation 
is complicated by the grain boundaries. We first sep- 
arate the defects into compact point defects and larger 
grain boundaries. For the point defects we determine 
the topological charge by taking the usual phase-angle 
path integral around the center of mass of the defect. 
Those defects with plus or minus unit charge are iden- 
tified as disclinations, while those with zero charge are 
dislocations. Then we can track the motion of each sin- 
gle defect. 



with Boyer and Vihals' discussion [7]. 

The most important motion of grain boundaries is that 
they can move over long distances and combine with 
other grain boundaries. As shown in Fig. 7, two grain 
boundaries can combine to form a larger grain bound- 
ary. Thus the number of grain boundaries decreases while 
their average size increases. This process happens on a 
time scale of the order 1000 dimcnsionless time units. 
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FIG. 6. A typical example of the track of a dislocation's "mass 
center", see appendix A. The dislocation moves along the arrow. 
It starts at t = 1590 and disappears at t = 29220. Notice the small 
arcs along the curve, the diameters of the arcs are about 2n, which 
is equal to two layer spacing. 

As an example of the method we show in Fig. 6 the 
path of a dislocation. We see that some dislocations 
travel over long distances during very long times. It 
seems that the dislocations are more stable when com- 
pared to grain boundaries and disclinations. Our simula- 
tions on 256 x 256 system show that after the annihilation 
of point defects and grain boundaries, some dislocations 
still exist in the system. Our simulations show that there 
are also dislocations which are pinned and move little. 

The number of disclinatins is quite small. And we no- 
tice that they are rather immobile, which is consistent 



FIG. 7. The combination of two grain boundaries in a 512 X 512 
system. The portion shown is 200 X 200. From left to right and top 
to bottom, the times are t = 2880, 2955, 3030, 3105, 3180, 4710. 
Not all the points in the grain boundaries are shown. 
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FIG. 8. The motion of a grain boundary in a 512 X 512 sys- 
tem. The portion shown is 200 X 200. From left to right and top 
to bottom, the times are t = 11415, 13665, 15915, 18165, 20415, 
22665, 24915, 27165, 29415. Again not all the points in the grain 
boundary are shown. 



As shown in Fig. 8, one grain boundary can sweep 
across a quite large area. At the same time its size de- 
creases. This process occurs on a time scale of the order 
10000. According to our observations, the grain bound- 
aries' motions also relieve the stripe curvatures through 
disclination annihilations. After one grain boundary 
passes through a disclination, the disclination disappears. 
This is consistent with Boyer and Vihals' prediction [7]. 
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FIG. 10. The average number of grain boundaries and the av- 
erage number of points in a grain boundary. 256 X 256 system 
averaged over 40 trials. 
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FIG. 9. The number of points in point defects and grain bound- 
aries in the 256 X 256 system. The data for all defects and domain 
walls are averaged over 40 trials. The others are averaged over 38 
trials. 



Next we focus on the statistics of the defects gener- 
ated by the model. In Fig. 9 we plot the total number of 
points in grain boundaries, dislocations and disclinations 
separately for the 256 x 256 system. We see that the grain 
boundaries dominate. In the scaling regime (t s < t < t c ), 
we see that the number of points corresponding to grain 
boundaries and all defect points, the curves a and b can 
be fit to ~ t -1 / 3 . At late stages the disclinations disap- 
pear, while the dislocations and grain boundaries persist. 
The number of disclinations decreases much faster than 
the other defects. 



In Fig. 10 we plot the average number of grain bound- 
aries n and the average size of a grain boundary I for the 
256 x 256 system. We use the number of points in one 
grain boundary as a measure of its size. For t s < t < t c , 
h decreases but I increases due to the combining of grain 
boundaries. The shrinkage of their sizes is not as im- 
portant as the combinations. However for t c ~ 9000 the 
correlation length L n is the same order as the system's 
size, and the large grain boundaries stop growing. After 
that the shrinkage is important [16] . In the scaling regime 
{t s < t < t c ), n ~ t~ 0A5 , and f~ t 013 . So hi ~ i" 1 / 3 , 
which is consistent with Fig. 9. 



6000 




'() 20000 ' 40000 

Time 



FIG. 11. The average number of points in defects and grain 
boundaries in the 512 X 512 system. Most of the points are in grain 
boundaries. All the curves can be fit to t~ y . The data for defects 
and domain walls are averaged over 57 trials, and the other data 
are averaged over 20 trials. 
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FIG. 12. n and I v.s. time in the 512 X 512 system. The data 
are averaged over 57 independent trials. 



In Fig. 11 we plot the total number of points in grain 
boundaries, dislocations and disclinations separately for 
the 512 x 512 system. The number of points in discli- 
nations is much smaller than that of dislocations and 
at very late stages it decreases to the order of 1, which 
in fact indicates the disappearance of disclinations. Now 
the scaling regime extends to much longer times. The do- 
main walls and dislocations' scaling exponents are both 
1/3, which is same to the scaling of the energy. How- 
ever, the scaling exponent of disclinations, given by 0.57 
is much larger. So we conclude that disclinations are not 
the dominant structures in SH system. 

In Fig. 12 we plot the number of grain boundaries n 
and the average size of a grain boundary I for 512 x 512 
system. All of our data falls in the scaling regime.. The 
plot of the average number h of grain boundaries ver- 
sus time t, can be fit to n ~ t~ QA9 and the average 
size I v.s. time t, can be fit by I ~ t° . So we have 
nl ~ L^ 1 ~ i -1 / 3 , which is consistent with the results 
shown in Fig. 11. 

Although we did not count the number of disclinations 
directly, it is proportional to the number of lattice points 
in disclinations, i.e. t~ ' 5 '. This is because the average 
number of lattice points in one point disclination, which 
is about 10 ~ 20, is quite stable during the simulation. 
By the same reasoning, we find that the number of dislo- 
cations is proportional to i™ 1 / 3 . It is interesting to note 
that the number of grain boundaries scales as t~ 0A9 . This 
exponent is near to that for disclinations. 

The number of grain boundaries is about 5 at t ~ 
70000, the number of disclinations is about or 1 at 
t ~ 50000, and the number of dislocations is on the order 
of 10 at t ~ 50000, as can be seen in Figs. 11 and 12. 



FIG. 13. A typical example of probability density P v.s. speed 
v for the speed distribution of dislocations at time 1350 for the 
512 x 512 system. The distributions at other times have approxi- 
mately the same shape. Averaged over 56 trials. 



Since we can track the motion of each defect, we can 
measure their speeds. We define the speed of each as the 
speed of its mass center (see appendix A for the definition 
of "mass center"). If in a time At, the mass center trav- 
els over a distance Ad, then the speed is v = Ad/ At. If 
At is small enough, we found that v = has the biggest 
probability. If At is large enough, all the details are 
coarse-grained and we observe a continuous distribution 
of the speed and the largest probability appears at a non- 
zero speed, as is shown in Fig. 13 where At — 60. We 
measured the speed distributions of domain walls and 
point defects separately. As we have already seen, for 
point defects the number of disclinations is much smaller 
than that of dislocations, so what we measure in the lat- 
ter case is in fact the speed distribution of dislocations. 



The speed distribution has a long tail which decreases 
as a power law. The numerical fits at different times give 
us different exponents. However the tail exponents at dif- 
ferent times do distribute in a narrow region, as is shown 
in Fig. 14. The exponents of grain boundaries are quite 
different from those of point defects. If we ignore the ex- 
ponents at very early times when grain boundaries just 
begin to form and at the very late times when the grain 
boundaries have already disappeared, the mean value of 
the grain boundaries' exponents is —1.50, and that of 
point defects' exponents is —2.10. 
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FIG. 14. The power law exponents of the speed distribution's 
tail at different times. Ignoring the data points at very early times 
and very late times, the mean value of the exponents is —1.5 for 
grain boundaries, —2.1 for point defects and —1.7 for all the defects. 
Averaged over 61 independent trials. 
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FIG. 15. The average speed for defects and grain boundaries. 
From top to bottom, the curves have the form vo X t~ x with vo be- 
ing a constant. From top to bottom, x = 0.35±0.01, x = 0.41±0.01 
and x = 0.48 ± 0.01. Averaged over 56 independent trials. 

We also measured the average speed of the point de- 
fects and grain boundaries as a function of time after the 
quench, as is shown in Fig. 15. The average speed of 
point defects decreases as ~ t~ QAS ; the average speed of 
grain boundaries goes as ~ i~ - 35 . The scattering of the 
points at late stages is due to the small data sample at 
those times.. 



VI. SUMMARY 



a control parameter of e = 0.1. We find in agreement with 
earlier workers that the kinetics in the ordering regime, 
before finite-size effects enter, arc dominated by the ex- 
istence of moving and coalescing grain boundaries. In 
this regime the average size of these grain boundaries is 
growing and they are relatively mobile. Under the influ- 
ence of finite size effects these grain boundaries shrink, 
the system becomes anisotropic and the ordering process 
speeds up. 

We also measured the speed distribution of all struc- 
tures that appear in the system. The average speed is 
decreasing as a power law and the distributions show a 
power-law behaviors at large speeds. However we can 
only get a rough estimate due to the poor statistics. 

Let us return to the question of whether the SH model 
gives a good description of the physical system studied 
by Harrison, et al. [2,3]. The SH model, for small control 
parameter e, does give coarsening with an exponent in 
roughly the same range as in the experiment (1/4 ~ 1/3). 
The ordering is constrained to be slower then the picture 
where one has a simple point defect pair annihilation 
process. However, the defect structures in the SH model 
and experiment appear quite different. The disclination 
quadrapolc annihilations seen in experiments are not ob- 
served in the late stages of the evolution of the SH model. 
In the SH model grain boundaries dominate the evolution 
in the scaling regime, but these structures appear to play 
a limited role in the experiments. We must make clear 
that our numerical results are for systems with many 
fewer roll periods compared to the experimental systems 
(10 2 compared to 10 5 ), so it is possible that things change 
as we increase the size of the ordering system. However, 
our study of 256 and 512 systems shows that they differ 
only in the time when the finite-size effect enters. This 
indicates that a even larger system will display the same 
behavior except that the finite-size effect enters at a even 
later time. So we conclude that the SH model does not 
give a physically faithful description of the ordering in the 
experimental system. This raises the provocative ques- 
tion: Are there many different types of scenarios for or- 
dering striped systems? We will address this question by 
looking at other competing models for striped formation 
elsewhere. 
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APPENDIX A: THE ALGORITHM FOR 
PICKING OUT DEFECTS AND FINDING 
DOMAIN WALLS 



We have studied the dynamics of the defect structures 
in the SH model after quenches to zero temperature with 



Hou, et al. [6] proposed a method, the HSG method, to 
measure the length of grain boundaries. They computed 
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the quantity A 2 = tp 2 + (Vi/>) 2 /qfi, and if the calculated 
A 2 is bigger than an upper filter 0.7 x Max(A 2 ) + 0.3 x 
AVG(A 2 ) or smaller than a lower filter 0.7 x Min(A 2 ) + 
0.3 x AVG(A 2 ), that point is counted as belonging to a 
domain wall. When e is small, this method gives quite 
good results. 

However, when e increases, the original filters are no 
longer applicable. They fail to pick out most of the 
points and the filters must be re-chosen. For example, 
at e = 0.6, the filters 0.5 x Max(A 2 ) + 0.5 x AVG(A 2 ) 
and 0.5 x Min(A 2 ) + 0.5 x AVG(A 2 ) can give a satisfying 
result; while for e = 0.75, OAx Max(A 2 )+0.6x AVG(A 2 ) 
and 0.4 x Min(A 2 ) + 0.6 x AVG(A 2 ) are the better 
choices. Sometimes this method is unable to pick out 
all the defects for any choice of filter. 

We introduce here a method which works for all e and 
picks out all of the defects and nothing more. 

First let us define some useful quantities. Suppose the 
system is discrete on the xy plane, with x = (square 
lattice). At a fixed time, starting from the order param- 
eter field ip(x) = V'iji we can define a director field n(x) 
as given by Eq.(ll) where V?/>(x) is defined by the usual 
finite difference scheme, i.e. 



of fields for our purposes. In the next step, we calculate 
</?(x) from sin<^(x) and cos</?(x) using 



(p(x) = arctan 



sin (p(x) 



cos tp(x) 



(A6) 



where we adopt the convention that — it < (p(x) < ir. 

In picking a filter we want to look at the spatial vari- 
ation of the v?( x ) field. V(/?(x) can be evaluated as for 
Vf/>(x). However, there is a subtlety here. For exam- 



ple, if <Pi-\-ij = tt — 5<j)i and (fi-i t 



5<j)2, where 



8(f>i and 8(j)2 are small angles, then the difference be- 
tween the nematic tensor Q a p(i + l,j) and Q a p(i — 
should be a small quantity. But (ipi + ij — <^j_i j j)/2 = 
tt — (5(f)i + 54>2)/2 ~ tt, which means that if we calcu- 
late the change rate of <p(x) in exactly the same way 
as Eq.(Al), we will get a wrong answer in this context. 
To avoid such a problem, we define the difference be- 
tween <p(x) and <p{x!) as the quantity with the small- 
est absolute value among the choices y(x') — <p{x) and 
^s(x') — ip(x) ± 2ir. And we use this quantity in deter- 
mining 



W(x) = 



ipi+i,j - 4>i-i,j ipi,. 



3+1 



^i,3 



2Ar 



2Ar 



(Al) 



where Ar is the lattice space of the system. In two- 
dimensional cases the nematic order parameter, Q a p is 
completely specified by the angle 



y>(x) = 2 0(x) 



where 



0(x) = arctan 



n x (x) 



(A2) 



(A3) 



Rather than using i/?(x) given by the two equations 
above we introduce some local smoothing. First we com- 
pute 



By = siny>(x) = 2n 2 .(x)n y (x) , 
B x = cos<p(x) = 2h x (x) 2 — 1 . 



(A4) 



Then we smooth these two fields using the iterative pro- 



cess: 



f(n+l)(i,j) = ^f(n)(i,j) + ^ 

(i'j')GNN 



/(„)(*',/) , (A5) 



where /(„) is simp or cos ip after n iterations, and JVJV 
means the 4 nearest neighbors of on the square lat- 
tice. This process will suppress the small fluctuations 
of y(x) away from the defects, while the variation of of 
<yj(x) near a defect core remains large. Our calculations 
show that 5 iterations provides a sufficiently smooth set 



A(x) = |Vc?(x) 



which is the key quantity in our analysis. 



(A7) 




ltU"lS* «.» " 



..................... ................. vl |, 

FIG. 16. In the lower graph, the vector field (cos ip(x), simp(x)) 
for the order parameter field shown above. The components of 
the vector field have been smoothed over 5 iterations. The lattice 
spacing is 7r/4, which means there are 8 points in one period of the 
layers. Not all the vectors on the lattice are shown. 
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Our method is based on the observation that at the 
core region of a defect (dislocation, disclination or part 
of a domain wall) the angle field y(x) changes rapidly, 
while in the region away from the defect's core the y(x) 
field is rather smooth, as can be seen in Fig. 16. Thus 
we can conclude with confidence that those points with 
larger change rates of y(x) must belong to some defect's 
core region or a part of a grain boundary. 




100 



FIG. 17. In the lower graph, the scalar field A(x) is plotted. 
This corresponds to the order parameter field shown above. A(x) 
is sharply peaked at the core regions ol the defects. 

We find that A(x) « away from defects, but increases 
very rapidly in the vicinity of any defect. An example 
is given in Fig. 17. Therefore as long as A(x) is large 
enough, we can identify the point x = (i, j) as part of the 
core region of a defect. Naturally we set up a threshold 
Aq, and any point with 4(Ar) 2 • A(x) > Ao is counted 
as belonging to some defect's core. Because the value 
of A(x) is much larger in the defects' cores than at any 
other places, a range of values of Aq can be used to find 
the positions of the defects core regions. With a smaller 
threshold the program will pick out more points in the 
core regions, and with a larger one it will pick out fewer 
points in the core regions. Our experience shows that if 
Ao takes the value of 2 ~ 10, the program picks out the 
same defects cores and grain boundaries. As is shown in 
Fig. 18 and Fig. 19, it picks out all the defects without 
irrelevant points. 




FIG. 18. Identification of all the defects in a 512 X 512 system 
(e = 0.1) with a threshold Ao = 3.5. At each defect core region, 
the A field for many points exceeds the threshold. The red points 
belong to domain walls, the green ones belong to dislocations and 
the blue ones belong to disclinations. 




FIG. 19. Identification of all the defects in a 512 X 512 system 
(e = 0.5) again with a threshold Ao = 3.5. Apparently the defect's 
density in this system is greater than the density in Fig. 18. The 
domain walls are much smaller for the system with larger e. There 
are no domain walls for e > 0.6. 

After we have used the above algorithm to pick out the 
points in the core regions of the point defects and grain 
boundaries, we can distinguish between these two struc- 
tures. The difference between them is obvious. The point 
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defects are compact in space while the grain boundaries 
are ramified. 

First, we must group the points we have identified ac- 
cording to whether or not they are in the same structure. 
The points in one point defect core or grain boundary are 
picked out because the director field changes drastically 
on those sites. They are very near to each other. How- 
ever, they may not be neighbors. So we define a filter do, 
and when any two points' distance is less than ao, they 
are supposed to be in the same defect or grain boundary's 
core region. We use the cluster multiple labeling method 
of Hoshcn and Kopclman [17] to pick out such point clus- 
ters. Thus given the system's status at any time, we can 
find those sets of points corresponding to each individual 
defect or grain boundary. 

Now we measure the approximate size of these struc- 
tures and then distinguish between point defects and 
grain boundaries. We use the number of the points in the 
set as the size of the corresponding structure. This ap- 
proximation reflects the actual size of the defect or grain 
boundary quite well. Then we define a filter l 0l and when 
the structure's size is larger than l , we regard the cor- 
responding structure as a grain boundary, otherwise it's 
taken to be a point defect (dislocation or disclination) . 

We employed ao — 5 Ar where Ar is the lattice spacing 
and lo = 18. The results are quite satisfying. 

After we have picked out the point defects, we can de- 
vide them into disclinations and dislocations. We follow 
Harrison's method [2,3]. Given the angle field <p(x) com- 
puted in Eq.(A6) after the smoothing process, we do an 
integral of the variation of </?(x) over a counterclockwise 
close path around the "mass-center" of a point defect, 
which is defined below. The condition for a defect to be 
a disclination is 

j^ds = ±2tt (A8) 

The integral is zero if the defect is a dislocation. To make 
the computation easier, we choose a 16 x 16 square with 
the mass-center at its center as the integration route. 

To record the motion of one single defect or grain 
boundary, we track the motion of the corresponding 
point set's "mass center", which is defined as follows. 
Suppose the point set has n points with coordinates 
Ti, i = 1, 2, n. Then the "mass center" of the point set 
is defined as r = ^r^/n, just like the usual mass cen- 

i 

ter definition in classical mechanics but with all masses 
equal to one. 

In the evolution of the SH model, we sample the sys- 
tem every 500 time steps (in our case this is equal to 
15 dimensionless time units), which is a quite short- 
time period in the simulation. We then identify all the 
dislocations, disclinations and grain boundaries, distin- 
guish among them, and compute their centers of mass. 
Suppose at time t\, we have the set of mass centers 



P = {pi, i = 1, 2, ...,m}, and at time t 2 , the "mass cen- 
ter" set is Q — {qj, j = l,2,...,n 2 }; usually m ^ n 2 . 
Define dpQ(i 7 j) = \pi — qj\. We assume that the defects 
and grain boundaries do not move much in such a short 
time period. So if there exist two integers k G [l,ni] and 
I G [1, n 2 ], such that 

d PQ (k,l)= min d PQ (k,j)= min d PQ (i,l) , (A9) 

je[l,n 2 ] i€[l,ni] 

it is quite reasonable to believe that and q; are just 
the same defect's or grain boundary's "mass center" at 
two successive times. Using this method, we are able to 
find out the trajectories of the "mass centers" as time 
goes on. Not all points in P and Q can be grouped into 
such pairs. On the one hand, this is because n\ ^ n 2 ; 
on the other hand, this is also due to the criterion (A9) 
applied onto p& and q;. Physically, this is consistent 
with the phenomena of the defect annihilation and the 
combination, split and shrinkage of grain boundaries. 



APPENDIX B: MEASUREMENT OF THE 
NEMATIC CORRELATION FUNCTION 

To probe the stripes' increasingly orientational order, 
we define the correlation function which is similar to the 
one employed by Christensen and Bray [10]. 

C nn (r, t) = cos [y?(x + r,t)- <p(x, t)\ ) 

X 

= ^2(cos(p(x + r,t) ■ cos <p(x,i)} + 

X 

+ ^H< sin ^( x + r > t )- sm ( x >*)> ' ( B1 ) 

X 

where N 2 is the area of the system and the angular brack- 
ets denote the statistical average over different initial con- 
ditions. The definition of the angle <p(x) is given in ap- 
pendix A. 

Now in Eq.(Bl) the function has been split into two 
parts which have the same form 

G ( r ' t ) = ^E(/( x+r ' t )'^')> - ( B2 ) 

X 

with /(x, i) = cosip(x,t) and f(x,t) = sin<^(x, t) sepa- 
rately. Eq.(B2) can be easily calculated by fast Fourier 
transformation (FFT). First FFT f(r,t) to obtain its 
Fourier components /(k,i). Then G(k,i) = ( |/(k,t)| 2 }, 
and inverse Fourier transformation gives G(r,t). We 
compute the two parts in Eq.(B2) separately, and then 
add to obtain the correlation function C nn (r, t). 
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